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Abstract 

Polymerization of proteins is a bioehimieal proeess involved in different diseases. Mathemati- 
cally, it is generally modeled by aggregation-fragmentation-type equations. In this paper we con- 
sider a general polymerization model and propose a high-order numerical scheme to investigate the 
behavior of the solution. An important property of the equation is the mass conservation. The 
fifth-order WENO scheme is built to preserve the total mass of proteins along time. 
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1 Introduction 

The central mechanism of amyloid diseases is the polymerization of proteins : PrP in Prion diseases, 
APP in Alzheimer, Htt in Huntington. The abnormal form of these proteins is pathogenic and has 
the ability to polymerize into fibrils. In order to well understand this process, investigation of the 
size repartition of polymers is a crucial point. To this end, we discuss in this paper the mathematical 
modeling of these polymerization processes and we propose numerical methods to investigate the 
mathematical features of the models. 

Mathematical models are already widely used to study the polymerization mechanism of Prion 
diseases [IIIIIIHIISTIIIHIIMIESIIMIEZIE], Alzheimer [HI [381 HH] or Huntington ^. Such models 
are also used for other biological polymerization processes [H [3] and even for cell division [21 [131 SS] 
or in neurosciences |47j . 

Another field where we find aggregation-fragmentation equations is the physics of aggregates (aerosol 
and raindrop formation, smoke, sprays...). Among these models (see [3l] for a review), one can mention 
the Smoluchowsky coagulation equation [201 [25l [261 [351 [H] with fragmentation [2T1 [22l [321 [331 [3T] 
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and the Lifshitz-Slyosov system [TlEKinilMllSlllllHSllle]. In[8l[29]a Smoluchowsky coagulation 
term is added to the Lifshitz-Slyosov equation. 

In this paper we are interested in a model including polymerization, coagulation and fragmentation 
phenomena. We consider a medium where there are monomers (normal proteins for instance) charac- 
terized by the concentration V{t) at time t and polymers (aggregates of abnormal proteins) of size x 
with the concentration u{x,t). The dynamics of the density function u{x,t) is driven by the system 

-V{t) = - T{V{t),x)u{x,t)dx, 

^u{x,t) = -^(T{V{t),x)u{x,t)^ +Q{u){x,t), (1) 

^it(0,t) = 0, u{x,0) = uo{x)>0 and ^(0) = Vb > 0. 

The monomers are attached by polymers of size x with the polymerization rate kon{x). Depolymer- 
ization occurs when monomers detach from polymers with a rate A:ofr(a^)- Hence the transport term 
writes 

T{V,x) = Vkon{x)-Ks{x). (2) 

The two functions /con and kog are piecewise derivable but can be discontinuous. They are positive 
except that kon can vanish at zero. In this case, or more generally when T{V{t),0) < 0, the boundary 
condition on u{0, t) is not necessary since the characteristic curves outgo from the domain. The choice 
of the bounbary condition u{0, t) = is justified later. 

The coalescence of two polymers and the fragmentation of a polymer into two smaller ones are taken 
into account by the operator Q. More precisely, denoting by an aggregate of size x we have 

kc{x,y) 

Ax + Ay > Ax+y coagulation 

, k{{x,y) A A 

Ax+y > Ax + Ay fragmentation. 

Thus the coagulation-fragmentation operator is Q = Qc — Qf with 



2 rx poo 

Qc{u){x) = - kc{y,x-y)u{y)u{x-y)dy-u{x) k^x , y) u{y) dy , (3) 
'0 Jo 
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Q{{u){x) = -u{x) kf{y,x-y)dy- k{{x , y) u{x + y) dy . (4) 
^ Jo Jo 

The coalescence of two polymers of size x and y occurs with the symmetric rate kc{x,y) = kc{y,x). 
This rate is a nonnegative function as the fragmentation symmetric rate fcf(x, y) = k{{y, x) with which 
a polymer of size x + y produces two fragments of size x and y. 



There is a difference between V{t) and u(x = 0,t). In biochemical polymerization processes, small 
polymers are very unstable and thus do not exist. When they appear by detachment from a longer 
polymer, they are immediately degraded into monomers. Thus, the quantity of small polymers vanishes 
while the quantity of monomers is very high. To reflect this in the mathematical model, a quantity 
V{t) of monomers is introduced, which is different from the quantity of small polymers u{x = 0,t). 
The evolution of the first one is given by an ODE while the second one is forced to be equal to zero 
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through the boundary condition u{0, t) = 0. A consequence of this distinction is that starting from 
Uo{x) = and Vq > there is no evolution : the concentration of monomers is constant in time, 
V{t) = Vq, and the concentration of polymers remains null, u{x,t) = 0. This is a very intuitive and 
natural behaviour which is important to preserve for biological applications. 

In the modeling, the distinction between V and u{x = 0) induces a separation of the polymerization- 
depolymerization process from the coagulation- fragmentation. Indeed the aggregation of a monomer 
to a polymer can be seen as a coagulation but the resulting polymer has same size x than the initial 
one, since a monomer is very small compared to the typical size of a polymer. So a transport term is 
more accurate to model this phenomenon than an integral term (see |14j). 

There is also the fact that when a small polymer is degraded into monomers, it increases the quantity 
of monomers. In a discrete model, this term appears in the equation on V (see ng in (4^). In 
the continuous model ([1]) this term can be neglected since the quantity of monomers produced by 
degradation of small polymers is very small compared to the total quantity of monomers. 



2 Mass Conservation 

The mechanism of polymerization is nothing but a rearrangement of the proteins, there is no creation 
and no disparition. So the total quantity of proteins has to be constant in time and this is the case in 
model ([1]). We define the total mass of the system as 



since a polymer of size x "contains x monomers". Integrating the equation on u{x,t) multiplied by x 
and adding the equation on V we obtain 



so the total mass is conserved along time. This is a very important property that we want to keep in 
the numerical scheme and for this we rewrite equation ([1]) under a conservative form. 

2.1 Conservative formulation 

The classical discretization methods for transport equations are mass preserving. So the idea is to 
write the coagulation- fragmentation operator Q, which preserves the mass, under a conservative form 
in order to use a transport scheme. For this we follow the paper [23] where such a transformation is 
made : 




(5) 



Vt > 




(6) 




where the operator C{u) is given by 




V 



(7) 
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and T{u) is 

px poo 

T{u){x) := / yk{{y,z)u{y + z)dzdy. (8) 

Jo Jx-y 

Under this form, the mass conservation is clearer and the use of conservative schemes possible. 

A useful consequence of the property ([6]) for the numerical scheme is that the ODE on V can be 
replaced by a mass conservation equation (see [29] ) 

poo 

Vt>0, V{t) = Vo+ j x{uo{x) -u{x,t))dx. (9) 

Numerically, this equation is much easier to compute than the ODE to be solved. Moreover ([9]) 
provides an explicit expression for V as a, function of u. So we set 

Q{u){x) := (Vq + y[uo{y) - u{y)\dy^kon{x) - Ka{x) (10) 

and we obtain a new equation equivalent to ([I]) 

x^u{x,t) + ^I^M^(x,i) + ^^ix,i) - ^^(^'*) = Giu)u{x,t), ^^^^ 

u{0,t) = 0, u{x,0) = uq{x). 

In this equation (jlip . we have written the transport term as 

d\g(u)u\ , ^ d\Q{u)xu\ , ^ ^, ^ , ^ , ^ 

X ^ U x,t) = ^ U x,t) - g{u)u{x,t). (12) 

This formulation enhances the relation 

d f°° f°° d 

— / xu{x,t)dx= / Q{u)u{x,t) dx = ——V{t) (13) 
dt Jq Jq dt 

and allows to preserve this property numerically when using conservative transport schemes. 
2.2 Domain truncation 

Numerically, equation (jll|) is solved on a truncated domain [0, B\ so the integration bounds have to 
be changed in order to keep the mass preservation. For the coagulation term, we introduce as in |23j 

f-x r-R-y 

C^{u){x) = 11 ykc{y,z)u{y)u{z)dzdy 

Jo Jx-y 
t-x rR 

= 11 ykc{y,z - y)u{y)u{z - y)dzdy, 

Jo Jx 

and for the fragmentation 



rx rR—y 

F^{u){x) := ykf{y,z)u{y + z)dzdy 

Jo Jx-y 
<-x i-R 

/ yhiy,z - y)u{z)dzdy. 

Jx 
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With this truncation, we have C^{u){0) = C^{u){R) = T^{u){0) = 7^{u){R) = 0. So the total mass 
does neither increase nor decrease with respect to time if we consider the coagulation and fragmentation 
processes. If we look at the effects of this truncation on the original coagulation and fragmentation 
operators we have 



1„ .T,, . 1 



R-x 



Qc{u){x):= dxC (u)(x) = - / kc{y,x -y)u{u)u{x -y)dy -u{x) kc{x , y)u{y) dy 

and 

Q^{u){x) := dxT^ {u){x) = -u{x) j kf{y,x-y)dy- h{x,y - x)u{y) dy. 

X Jo Jx 

In the coagulation term, the truncation corresponds to the assumption that a polymer of size x cannot 
coagulate with a polymer of size greater than R — x. Concerning the fragmentation term, it is nothing 
but the assumption that polymers of size greater than R cannot split. Biologically they are the natural 
assumptions to avoid the loss of mass. 
Concerning the transport term, the only way to avoid the loss of mass is to set 

g^{u){R,t) =0. (14) 

The meaning we give to this relation in the numerical scheme is exposed in Section 13.21 It is useless 
to do such a truncation for x = since xu{x, t) vanishes when x = 0. 
Finally we obtain a conservative truncated equation for x G (0, i?) 

t) + "iS^f""] t) + ^^^(X, t) - ^^^(X, t) = 5(.h)u„(x, t), ^^^^ 

^UR{0,t)=0, ur{x,0) = uo{x). 

When there is no transport term, convergence of the solution of Equation (jl5p to the solution of 
Equation ([1]) when i? — ?> oo is proved in \17\ \3^ [32} \3T\ [57] under growth conditions on kc and k{. 



3 A Fifth Order WENO Scheme 

In order to obtain a mass preserving scheme, we consider equation (jlip as a transport equation 
and for high order accuracy we choose a fifth-order WENO (Weighted Essentially Non Oscillatory) 
reconstruction for the fluxes. This high order scheme is comonly used |12[ since it is not more 
complicated to implement than a third order WENO one for instance. 



3.1 Numerical fluxes 

Before using the WENO reconstruction we have to know if the fluxes are positive or negative in order 
to appropriately upwind the scheme. Concerning the coagulation and the fragmentation terms, we 
consider a positive upwinding as suggested in [23]. For the transport term [t7(n)xti] we have to 
make a flux splitting because Q has no sign. A natural splitting here is to put the terms of G that are 
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preceded by a plus sign in the positive part and the terms preceded by a minus sign in the negative 
part, namely Q = + where 



gf{u){x) = [Vq + yuo{y)dyjkon{x), ^^^^ 
Giiu){x) = -(^J^ yu{y)dy^kon{x) - kos{x). 



An other decomposition is the polymerization-depolymerization one 

Soiu)ix) = (Vo + J^y[uoiy) - u{y)]dy^konix), 
gQ{u){x) = -kos{x). 



(17) 



The term g^ is necessarily positive because Vo + /o°° y [uo{y) - u{y, t)] dy = V{t) > 0. With these two 
flux splittings, we built others by convex combination. For any A € [0, 1] we set gx = At/i + (1 — A)t/o 
which gives 



g^{u){x) = {Vo + y[uo{y) - u{y)]dy + X J^^ yu{y)dyjkon{x), ^^^^ 
Sxiu)ix) = -X(j^ yu{y)dy]kon{x) -kosix). 



We also consider the Lax-Friedrichs scheme which corresponds to 

r gMu) = ^{g{u) + m), 

with m = maXa;>o |^('")|- This term has to be computed at each time step because g{u) depends on 
time. 

Finally, the WENO reconstruction is done with the fluxes 

H+{u) = g+{u)xu + C{u) - F{u), 

(20) 

H (u) = g {u)xu, 



and the choice between the different flux splittings is discussed in Section 14.21 
3.2 WENO reconstruction 

The point of view adopted here is the finite difference one, as recommanded in [55], because it is 
better than the finite volume in terms of operation counts. We assume the spatial domain [0, R] is 
divided into uniform cells and we denote Xi = iAx for < i < with Ax = ^. We use the 
WENO formulation of Jiang and Peng [53] which consists in applying WENO to approach the spacial 
derivative directly on the nodes of the grid. The spatial derivative dx{H~^ (v) + H" (v)) is approximated 
at the point Xi by 

1 



Ax 



H. I + , — H. ^ — H. 1 



where the fifth order accurate numerical flux H'^ i is given by the WENO reconstruction. For each 

node Xi we denote by the numerical approximation of H^{v{xi)). The stencil choice for each flux 
is specifled in Figure [Tl and the fluxes H^^ ^ are expressed as convex combination of the Hf^ of the 

stencil. Let us detail how we proceed : 
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for H- J 
for H7 J 
for H+ J 
for J 
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Figure 1: stencil choice 
For the regularity coefficients we define for each previous flux 



Si = 


'Aw, 


- 2W2 + W3 


' + \{W,- 


AW2 + 31^3)^ 


S2 = 


12^ 


- 2W3 + W4 






S3 = 


12^ 


- 2W4 + W5 


2 + i(3VF3 


- 4W4 + W^f 


hts 











Wr 



dr ,3,6,1 

with Qr = Ml = , do = , d-i = 

e + Sr 10' 10' 10 



where e is introduced to prevent the denominator from vanishing. 
Finally we take the different flux parts given by 



H 



W3 7W2 llWA 



1,2,3 



. , , , -W2 5W, WA , fWi , 5Wi W5 



Wl 



Wi 1W2 iiWsA (-^2 WA (W3 5W4 

-^^ — H — ^ I + ^^2 ( — — V + ^ + ^3 + 
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6 J 



6 



6 
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6 
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Concerning the boundaries x = and x = R,we compute the fluxes using the WENO reconstruction 
with ghost points x_3, X-2, X-i, and xn+i, xiy^2, ^^tv+s. In the first three points we use the boundary 
condition n(0, = to set u_3 = ii_2 = u_i = and for the last three ones we use G{u){x = R) = 

to put g^,. = , 2 = g^,^ = 0. 
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3.3 Integration method 



For the integral terms, we use a fifth order composite rule already used in |54] . If denotes an 
approximation of f{xk), the method can be written as 



k=i 



where 



2^ Ik - + 240 ^'^^ 240 720 "^^^^ 



k=i 



739 211 299 251 



j-4 



fe=i+4 



if j — i > 6. This method is based on polynomial interpolations of the function /. 

On the first interval, we integrate without using the boundary value /o = because the solution 
can be discontinuous at x = 0. So we use the fifth accurate approximation 

V/. _ 55 59 37 9 

k=0 

Finally for the intervals at the boundaries we have 

^ ' f': = g/i - + g/s - 2/4, 





3 



^, 21 9 15 3 

Z_j Ik - y Ji - g/2 + + g J4 - ^/5, 

Ik - -^II - g/2 + Y-'3 + + Ij^h - ^/6, 
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Ik — -^II - ^/2 + -^/3 + -^H + + - ■^/7i 



Vf -^^f ^ ^2 25 1 1 

Ik — -^II - ^J2 + -^/3 + + /5 + + - ^/S) 



g.^ g , , , 24-'" 2-" 24" 

N y 

2^..v 24-'^' 24" 



>p , 9 19 5 1 



iV-1 
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' /fc = g/iV + 3/^'-! + 

N-2 

^ 1 31 7 13 1 

' /fc = -^In + + 8-^^-2 + - 24 

iV-3 

N-4 

^ 1 31 5 25 25 1 1 

We use this quadrature method to discretize the operators and J-^ with 

i AT 

Cf - J-f = (Ax)2 J] ' ' (A;,^,z_,n,n;_,. - k^i.^ui) (21) 

j=0 «=i+l 

as suggested in [23] . Grouping the two terms in an unique summation is lighter regarding to operation 
counts. 

3.4 Time discretization 

The time step is denoted by At and changes along time because of the CFL stability condition that 
is time dependent. For the time discretization we choose a third order Runge-Kutta method. To 
approach the time evolution of an equation dtu = L{u), we compute at time nAt 

=u" + AtL(n") 

and 

^(2) = 3^n^l^(i) + lAtL(n«), 
4 4 4 

where is an approximation of u{nAt). Then the approximation of v at time (n + l)At is given by 



3 3 3 ^ ' 



This method is an explicit one, so to ensure the stability we compute the time step At at each iteration 
thanks to the CFL condition 

At < min{(G + C + F)^i} (22) 

where 

G = ^ sup (e+ -Gr), <^ = sup I ^ ' kljUj I and ^ = sup | ^ ' k^^^j | • 
For instance the Lax-Friedrichs decomposition leads to Glf = m/Ax. 
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4 Numerical Simulations 



4.1 Parameters 

Numerical values for the polymerization and fragmentation rates can be found in the biological liter- 
ature (see [Ml [521 SO] for instance). It is of importance to note that the models considered in these 
papers are discrete ones, so some computations are necessary to deduce adimensional numerical values 
for the continuous model ([T|). Another point is that the parameters of these models do not depend on 
the size of polymers, so we can only obtain mean values for the size-dependent parameters. 

We choose to use the values of the recent paper [30] to do numerical simulations. The mean length of 
polymers for the initial distribution is estimated to be 1380. With the continuous model we reduce this 
value to 0.2 by considering an initial profile equal to a positive constant on [0, 0.4] and null for x > 0.4 
(see the first plot of Figure [5]). Thus we define a parameter e = 0.2/1380 ~ 1,4 x 10""^ which allows 
to go from a discrete model to a continuous one (see |14j for more details). The values we find are 
for instance 2.9 x 10~^^M~^s~^ for the polymerization rate and 2.1 x lO^^s^^ for the fragmentation 
(where M represents the concentration in mol and s the time in second). The polymerization rate 
appears in a derivative term, so the value of the discrete model has to be multiply by e to obtain the 
continuous accurate value 4 x 10~^;uM~^s~^. Conversely, the fragmentation rate which appears in an 
integral term has to be divided by e and we find 1.5 x 10~^s~^. Concerning the depolymerization and 
coagulation, they are neglected in the models of [39^ [52l I30j . So we consider numerical values that 
seem to be reasonable compared to the previous ones. 

In the present study, the parameters are assumed to be size dependent as suggested in [U [6] and 
their choice is now presented and motivated. Concerning the numerical coefficients, they are chosen 
in order to have mean values of the same order than the values previously obtained from |30j . 
For the polymerization we assume that small polymers have a different behavior compared to the big 
ones. We consider a critical size Xc = 0.5 such that polymers of size x < Xc convert monomers with 
the rate 

A:W(x) = (4X + 0.2) x 10~^nM~^s"^, 
and for x > x^ with a constant rate 

A;(2)(x) = 4 X IQ-^^xM-^s-^. 

For the fragmentation kernel, we use the classical assumption that the fragmentation probability 
depends only on the size x + y of the polymer and we set 

10 + (x + y) 

The depolymerization is assumed to be constant and of the same order than the fragmentation. We 
discuss the dependence on this parameter by considering different intensities 

Ke{x) = r]X IQ-^s-^ with 2 < < 8. (23) 

Concerning the coagulation kernel, we do not use a classical one. Even if there is no space in model ([T]), 
we use a kernel which reflect some space effects. Indeed we consider that small polymers are very mobile 
and that the big ones, plaques, are very attractant. So the coagulation occurs preferentially between 
big and small aggregates. The kernel we choose is of the form 

, 3 



4|x — 2/1 2 
l + ix + y) 



kc{x,y) = X 10~*^^M~^s" 
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This kernel satisfies the growth assumption that we can find in [T7] to ensure the convergence of the 
solution when i? ^ oo if we consider the only coagulation-fragmentation process. 



coagulation (in jiM ^) 



fragmentation {in h ^) 





Figure 2: Profiles of the coagulation and fragmentation kernels. 

rR 

In [30] we also find numerical values for the initial datas Vq = 98 /xM and / xuo{x) dx = 0.21/xM. 

Jo 

This last value and the fact that the initial distribution of polymers is assumed to be under the form 
uq = est X ]l[o,o.4] lead to 

2.6 if < a; < 0.4 
if X > 0.4. 



uo{x) 



For the following simulations, the discretization is made on a domain [0, 5] with a number of nodes 
N = 200, so the mesh size is Ax = 0.025. 



4.2 Choose among the different flux sphttings 

First we deal with the CFL condition. Thanks to a triangular inequality, we obtain that Glf < Gq. 
Moreover, there is an explicit expression for Gx 



°5 = i;T{(''°+^ 



N 



N 

XjVL^ + (2A - 1) Ax ^ ' x,u 



off 



which shows that G\ increases with A. So we have that ifO<A<A<l then at each time step 
< G\ < G\. Notice also that, with the numerical values we have chosen, the quantity of polymers 



N 



= 1 ^j'^j 



increases with n. Indeed we can see on the Figure [7] that the quantity of monomers 



yn V{nAt) defined by the mass conservation V"' + Ax X^^Li' xju^ = Vq + Ax ' 



Xju'j decreases. 



The consequence on the CFL condition is that G^ increases with n if A > ^, decreases if A < ^ and 
is time independent when A = ^. Thus, regarding to the numerical computation, the fastest scheme 
is the Lax-Friedrichs one and then the computation time increases significantly with A. 



Let us now turn to the effects of the flux splitting on the size distribution. First we consider a 
depolymerization corresponding to = 8 in (|23|) and investigate the differences between the solutions 
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associated to the decompositions ^lfj Go and Gi. We can see on the Figure [3] that the solutions for 
Qo and Qi are close together for small times and then the behavior of Gq becomes closer to the Lax- 
Friedrichs one. The less oscillating scheme for t = 6h is the the Lax-Friedrichs one, but it is also the 
most oscillating at time t = 12h. Finally the solutions associated to the three flux decompositions are 
quite similar and they all present oscillations at some times, so we do not find with this simulation 
any reason to discard one of them. 



t = 6 t = 12 




0.2 0.4 0.6 0.8 1 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 



X X 
Figure 3: Comparison between the flux splittings Qi,f, Go and C/i for 77 = 8 x lO^^s 



If we change the depolymerization rate by considering kog = 6 x 10" s~ , we remark that the 
Lax-Friedrichs scheme is unstable (see Figure U]) while Go is stable. If we continue to decrease rj, we 
find with = 2 x 10~^s~^ that the ^o-scheme becomes unstable while G0.2 is stable. Thus the 
Lax-Friedrichs scheme and the ^^-schemes with A small has to be avoid to ensure stability when small 
depolymerization values are considered. 



o 





Figure 4: Unstability of some schemes when 77 decreases. Left : C/lf becomes unstable for 77 = 6. Right : Qq 
becomes unstable for 77 = 2. 
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Knowing that, we compare different stable schemes, namely Qx with 0.2 < A < 1. We can see 
on Figure [5] that for large times (t = 20h), the three flux decompositions provide a good behavior 
where there are strong variations of the solution. These locations are x = because of the boundary 
condition which enforces u{0, t) to vanish, and x = 0.5 because the transport term fcon is discontinuous 
at a; = 0.5. If we look at smaller times {h = 6h for instance) we can see that the larger A is, the less 
oscillating the curves are. But, as we already remarked, the quantity Gx is higher for A close to 1 and 
it increases with time when A > ^. Thus it is penalizing for the computation time to use high values 
of A. A good compromise could be to choose A = i since Gi does not depend on time. The other 

^ 2 

solution is to adapt the A when we change the parameters. 



t = 6 



t = 20 






120 - 




100 - 




80 - 








60 - 












40 




20 I 




- 




Figure 5: Comparison of the behavior of the solution for different A with 77 = 5. 



4.3 Interpretation of the numerical results 

We have considered that the mean size of the polymers at the initial time t = was 1380. This size 
can be multiplied by 5 along the polymerization process (see Figure [6] keeping in mind that the mean 
size is represented by 0.2 in this continuous model). So if we want to solve the discrete model, we 
have to consider a system of dimension close to 5000, and the computations are very heavy. If we 
limit this value to 200 keeping the discrete model, then we lose a lot of precision. It is the same 
for the continuous model if it is discretized with a first order scheme. That is why we use a fifth 
order discretization, and we can see the difference on Figure [6] : the fifth order scheme is able to 
capture strong variations of amplitude while the first order fiattens them. We also remark that the 
size repartition converges to a bimodal distribution, as observed by [56]. This asymptotic profile can 
be seen as an eigenvector of the operator Q — dxT (see [41^ I16j). 

The evolution of the quantity of monomers V{t) is plotted in Figure [7] for different values of the 
depolymerization rate kos- This quantity decreases since the monomers aggregate to polymers. Thus 
the mass of polymers increases and the speed of this evolution is similar to those observed by [30] . 
Concerning the dependence on k^g, the difference between the three curves is more significant when 
the time increases. For small times, when V{t) is close to Vq = 98/iM, the depolymerization can 
be neglected since is small compared to the product V{t)kon{x). Conversely, the equilibrium is 
reached when 4fV{t) = 0, so when k^s — Vko^ (see Equation ([T])). That is why variations of r/ infiuence 
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essentially the ratio between the quantity of monomers and the mass of polymers at the equilibrium 
as we can see on Figure [71 



5 Conclusion and future work 



We have written a high order conservative scheme for a polymerization-type equation. The choice 
of the flux splitting for the transport term has been discussed but the accurate decomposition remains 
unclear. Moreover the simulations show some oscillations for any choice of the flux decomposition. 
These points remain to be investigated for a better understanding and improvement. 

As we have remarked in Section 14.31 the size distribution converges toward an equilibrium which 
corresponds to an eigenvector. The fifth-order WENO scheme presented in this paper could be used to 
numerically compute such eigenvectors. Another application of the code is to solve inverse problems 
(see |15l I50j ) in order to determine the size dependence of the different parameters. 
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Figure 7: Evolution of tlic quantity of monomers for different depolymerization rates, with the scheme Qo.2- 
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